# Make LaTeX chunk for Table 1
install.packages('PRROC')
library("PRROC")
library("pROC")
# Load Models and Store AUCs
countries <- c('indo', 'colo')
dvs <- c("any", "high", "spike")
algos <- c('lasso', 'rf', 'gbm', 
           'nn', 'ebma')
metrics <- c('weight', 'acc_acc', 
                 'sens_acc', 'spec_acc',
                 'acc_sens', 'sens_sens',
                 'spec_sens', 'mse', 'prauc',
                 'auc', 'dv_mean')
performance <- array(NA,
                     dim = c(length(countries),
                           length(dvs),
                           length(algos),
                           length(metrics)),
                   dimnames = list(countries,
                                   dvs,
                                   algos,
                                   metrics))
for (country in countries) {
  for (algo in algos) {
    modelname = algo
    for (dv in dvs) {
      filename <- paste(modeldir,"/",
                    country,
                    "_",algo,"_",
                    dv,
                    "_full_mse.RData",
                    sep = "")
      load(filename)
      model <- get(paste(modelname,".results",sep=""))
      T <- length(model)
      metrics_by_year = array(NA,
                              dim=c(length(metrics), T),
                              dimnames = list(metrics, NULL))
      weights <- array(NA,
                       dim=c(length(algos[-length(algos)]), T),
                       dimnames = list(algos[-length(algos)], NULL))
      for (i in 1:T) {
          predictions <- as.vector(model[[i]]$fit.oos)
          predictions_acc <- predictions > model[[i]]$dt
          predictions_sens <- predictions > model[[i]]$dtsens
          if (algo == 'ebma') {
            predictions_acc <- predictions > model[[i]]$optim.dt
            predictions_sens <- predictions > model[[i]]$optim.dtsens
          }
          realizations <- as.vector(model[[i]]$actual.oos)
          ### EBMA Weights
          if (algo=='ebma') {
            weights[algos[-length(algos)], i] <- model[[i]]$optim.wts[paste0(algos[-length(algos)],
                                                                             ".weight")]
          }
          
          ### Performance Maximizing Accuracy
          metrics_by_year['acc_acc', i] <- error.rates(realizations, 
                                                   predictions_acc)['success.rate']
          metrics_by_year['sens_acc', i] <- error.rates(realizations, 
                                                       predictions_acc)['tp.rate']
          metrics_by_year['spec_acc', i] <- error.rates(realizations, 
                                                       predictions_acc)['tn.rate']
          
          ### Performance Maximizing Sensitivity
          metrics_by_year['acc_sens', i] <- error.rates(realizations, 
                                                       predictions_sens)['success.rate']
          metrics_by_year['sens_sens', i] <- error.rates(realizations, 
                                                        predictions_sens)['tp.rate']
          metrics_by_year['spec_sens', i] <- error.rates(realizations, 
                                                        predictions_sens)['tn.rate']
          
          ### MSE
          metrics_by_year['mse', i] <- mean((realizations - predictions)^2)
          
          ### PRAUC
          metrics_by_year['prauc', i] = pr.curve(scores.class0 = predictions[realizations==1],
                                               scores.class1 = predictions[realizations==0])$auc.integral
          
          ### AUC
          metrics_by_year['auc', i] <- roc(response = realizations,
                                           predictor = predictions)[['auc']]
          
          ### DV Mean
          if (algo == 'lasso') {
            metrics_by_year['dv_mean', i] <- mean(realizations) 
          }
      }
      performance[country, dv, algo, metrics] <- apply(metrics_by_year, 1, mean)
      if (algo == 'ebma') {
        performance[country, dv, algos[-length(algos)], "weight"] <- apply(weights, 1, mean)
      }
    }
  }
}

# Build Table with Labels
country.labels <- c(indo='\\textbf{Indonesia} (social conflict)',
                    colo='\\textbf{Colombia} (attacks and clashes)')
outcome.labels <- c(any="(a) Any violent event",
                    high="(b) $\\ge$ 5 violent events",
                    spike="(c) $\\ge$ 1 s.d. increase in events")
model.labels <- c(lasso='LASSO',
                  rf='Random\\\\Forest',
                  gbm='Adaptive\\\\Boosting',
                  nn='Neural\\\\Network',
                  ebma='EBMA')
metric.key <- c(weight = 'EBMA Weight', 
                   acc_acc = 'Accuracy', 
                   sens_acc = 'Sensitivity',
                   spec_acc = 'Specificity',
                   acc_sens = 'Accuracy', 
                   sens_sens = 'Sensitivity',
                   spec_sens = 'Specificity',
                   mse = 'MSE (Brier Score)',
                   prauc = 'P-R AUC',
                   auc = 'AUC',
                   dv_mean = 'Dep. Var. Mean'
)
breaks <- c('Threshold maximizes accuracy'=2,
            'Threshold maximizes sensitivity, while accuracy above 50\\%'=5,
            8)
table <- performance
dimnames(table)[[1]] <- country.labels[dimnames(table)[[1]]]
dimnames(table)[[2]] <- outcome.labels[dimnames(table)[[2]]]
dimnames(table)[[3]] <- model.labels[dimnames(table)[[3]]]

# Print to Latex

formatted_tab <- format(round(table, 3),
                        digit=3)
formatted_tab[formatted_tab == "   NA"] <- ""

printA1(table = formatted_tab,
        breaks=breaks,
        filepath = "tables",
        file = "table_A2",
        metric.key = metric.key,
        head.width = '1.2',
        head.height = '1.2')

